# ===============================================
# Paris ratchet paper
# Code for developing country growth rates
# ===============================================

#### Initialize 

#### Table of contents ####

#' 1. Load targets data
#' 2. Reshape data to calculate growth rates
#' 3. Get summaries


#### Initialize #### 

library(tidyverse)
library(readxl)


## Set working directory


# Define core ISO3C codes to use in analysis
core_iso3c <- c("AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEL", "BEN", "BFA", "BGD", "BGR", "BHR", "BHS", "BIH", "BLR", "BLZ", "BOL", "BRA", "BRB", "BRN", "BTN", "BWA", "CAF", "CAN", "CHE", "CHL", "CHN", "CIV", "CMR", "COD", "COG", "COK", "COL", "COM", "CPV", "CRI", "CUB", "CYP", "CZE", "DEU", "DJI", "DMA", "DNK", "DOM", "DZA", "ECU", "EGY", "ERI", "ESP", "EST", "ETH", "FIN", "FJI", "FRA", "FSM", "GAB", "GBR", "GEO", "GHA", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GTM", "GUY", "HND", "HRV", "HTI", "HUN", "IDN", "IND", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAO", "LBN", "LBR", "LBY", "LCA", "LIE", "LKA", "LSO", "LTU", "LUX", "LVA", "MAR", "MDA", "MDG", "MDV", "MEX", "MHL", "MKD", "MLI", "MLT", "MMR", "MNE", "MNG", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAM", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRU", "NZL", "OMN", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRK", "PRT", "PRY", "QAT", "ROU", "RUS", "RWA", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SOM", "SRB", "SSD", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SYC", "SYR", "TCD", "TGO", "THA", "TJK", "TKM", "TLS", "TON", "TTO", "TUN", "TUR", "TUV", "TZA", "UGA", "UKR", "URY", "USA", "UZB", "VCT", "VEN", "VNM", "VUT", "WSM", "YEM", "ZAF", "ZMB", "ZWE")

# Define a "max" function that works with NA
my.max <- function(x) ifelse( !all(is.na(x)), max(x, na.rm=T), NA)

## Load data

data_for_models <- readRDS("output/data_for_models20230920.Rds")


#### Load targets data from scratch ####

#' Most of this code is carried over from "00_make_analysis_data.R"
#' In the brace, loads targets and emissions data
#' Then calculates emissions growth rates by country

{
  # Load first NDCs
  ndc2015 <- read_csv("input/ndc_2015.csv") %>% 
    filter(iso3c %in% core_iso3c) %>% 
    arrange(iso3c) %>% 
    # Make proper date variables
    mutate(date = as_date("2015-12-12")) %>% 
    relocate(date, .after = ("country")) %>% 
    # Remove 'MtCO2' from values
    mutate_at(vars(c(cw_ndc1_ghgtarget_scenario_referencelevel,
                     cw_ndc1_ghgtarget_fixedlevel_goal,
                     cw_ndc1_ghgtarget_fixedlevel_goal_stretch)), ~parse_number(.)) %>% 
    mutate(cw_ndc1_ghgtarget_type = 
             ifelse(cw_ndc1_ghgtarget_type == "Base year", 
                    "Base Year", 
                    ifelse(cw_ndc1_ghgtarget_type == "Fixed level", 
                           "Fixed Level", cw_ndc1_ghgtarget_type)))
  
  # Load second NDCs
  ndc2021 <- read_excel("input/ndc_2021.xlsx", sheet = "2021 Dataset") %>% 
    arrange(iso3c) %>% 
    # Make proper date variables
    mutate(date = lubridate::as_date(submission_date_yy_mm_dd)) %>% 
    relocate(date, .after = ("country")) %>% 
    select(-submission_date_yy_mm_dd) %>% 
    # Remove 'MtCO2' from values
    mutate_at(vars(c(cw_ndc2_ghgtarget_scenario_referencelevel,
                     cw_ndc2_ghgtarget_fixedlevel_goal,
                     cw_ndc2_ghgtarget_fixedlevel_goal_stretch,
                     cw_ndc2_ghgtarget_fixedlevel_conditional)), ~parse_number(.))
  
  # Load GHG series
  primap <- read_csv("input/primap_v231.csv") %>% 
    filter(entity == "KYOTOGHG (AR4GWP100)") %>% # Aggregate using Kyoto GWP rules
    filter(`scenario (PRIMAP-hist)` == "HISTTP") %>% # Prioritize third-party emissions data over country data
    filter(`category (IPCC2006_PRIMAP)` == "M.0.EL") %>%  # Include all sectors
    pivot_longer(names_to = "year", values_to = "ghg_primap", `1980`:`2019`) %>% 
    select(iso3c = `area (ISO3)`, year, ghg_primap) %>% 
    filter(year > 1989) %>% 
    filter(!iso3c %in% c("ANNEXI", "AOSIS", "BASIC", "EARTH",
                         "EU28", "EU27BX", "NONANNEXI", "UMBRELLA")) %>% 
    filter(iso3c %in% core_iso3c) %>% 
    mutate(year = as.numeric(year)) %>% 
    mutate(ghg_primap = ghg_primap/1000) # into MtCO2
  
  # Get LULUCF status for countries' targets
  lulucf_flag <- read_csv("input/cw_ndc_data_highlevel_20220928.csv") %>% 
    select(iso3c = ISO, ndc_date, coverage_sectors_label) %>% 
    filter(iso3c %in% core_iso3c) %>% 
    group_by(iso3c) %>% 
    mutate(count = row_number()) %>% 
    filter(count == max(count)) %>%
    select(-count, -ndc_date) %>% 
    summarize(iso3c, 
              lulucf_flag = str_detect(coverage_sectors_label, 
                                       "including LULUCF"))
  
  # Load in WRI/CAIT/Climate Watch data
  cw <- read_csv("input/cw_ghg_20220928.csv") %>% 
    filter(Sector %in% c("Total excluding LUCF", "Total including LUCF"),
           Gas == "All GHG") %>% 
    pivot_longer(names_to = "year", values_to = "cw_ghg_", `1990`:`2019`) %>% 
    pivot_wider(names_from = "Sector", values_from = "cw_ghg_", ) %>% 
    select(iso3c = Country, year, 
           ghg_cw_ex  = `Total excluding LUCF`, 
           ghg_cw_lucf = `Total including LUCF`) %>% 
    mutate(year = as.numeric(year)) %>% 
    filter(iso3c %in% core_iso3c) %>% 
    as.data.frame() %>% 
    full_join(., lulucf_flag, by = 'iso3c') %>% 
    mutate(ghg_cw_chosen = ifelse(lulucf_flag == T, ghg_cw_lucf, ghg_cw_ex))
  
  # Merge GHG data and 2015 NDCs to get target levels
  targets2015 <- full_join(primap, cw, by = c("iso3c", "year")) %>% 
    full_join(., ndc2015 , by = "iso3c") %>% 
    arrange(iso3c, year) %>% 
    # Get emissions in base year and GHG 2010 as own variables
    mutate(ndc1_baseyear_emissions = ifelse((cw_ndc1_ghgtarget_type == "Base Year" &
                                               year == cw_ndc1_ghgtarget_referenceyear),
                                            ghg_cw_chosen, NA),
           ghg_cw_chosen_2000 = ifelse(year == 2000, ghg_cw_chosen, NA),
           ghg_cw_chosen_2010 = ifelse(year == 2010, ghg_cw_chosen, NA),
           ghg_cw_chosen_2015 = ifelse(year == 2015, ghg_cw_chosen, NA)) %>% 
    group_by(iso3c) %>% 
    mutate(ndc1_baseyear_emissions = my.max(ndc1_baseyear_emissions),
           ghg_cw_chosen_2000 = my.max(ghg_cw_chosen_2000),
           ghg_cw_chosen_2010 = my.max(ghg_cw_chosen_2010),
           ghg_cw_chosen_2015 = my.max(ghg_cw_chosen_2015)) %>% 
    ungroup() %>% 
    # Get absolute emissions levels in targets
    mutate(ndc1_absolute_uncond = 
             ifelse(cw_ndc1_ghgtarget_type == "Fixed Level", cw_ndc1_ghgtarget_fixedlevel_goal, # For fixed level targets
                    ifelse(cw_ndc1_ghgtarget_type == "Scenario", # For scenario targets
                           cw_ndc1_ghgtarget_scenario_referencelevel*(1-cw_ndc1_ghgtarget_scenario_goal/100),
                           ifelse(cw_ndc1_ghgtarget_type == "Base Year", # For base year targets
                                  ndc1_baseyear_emissions*(1-cw_ndc1_ghgtarget_baseyear_goal/100), NA)))) %>% 
    mutate(ndc1_absolute_cond = 
             ifelse(cw_ndc1_ghgtarget_type == "Fixed Level", cw_ndc1_ghgtarget_fixedlevel_conditional, # For fixed level targets
                    ifelse(cw_ndc1_ghgtarget_type == "Scenario", # For scenario targets
                           cw_ndc1_ghgtarget_scenario_referencelevel*(1-cw_ndc1_ghgtarget_scenario_goal_conditional/100),
                           ifelse(cw_ndc1_ghgtarget_type == "Base Year", # For base year targets
                                  ndc1_baseyear_emissions*(1-cw_ndc1_ghgtarget_baseyear_goal_conditional/100), NA)))) %>% 
    # Get emissions levels in targets as percentage of 2010 emissions (same benchmark as PRIMAP)
    mutate(ndc1_percentage_uncond = ndc1_absolute_uncond/ghg_cw_chosen_2010,
           ndc1_percentage_cond = ndc1_absolute_cond/ghg_cw_chosen_2010) %>% 
    # Get an average across both NDC targets (unconditional and conditional)
    rowwise() %>% 
    mutate(ndc1_percentage_average = mean(c(
      ndc1_percentage_uncond, ndc1_percentage_cond), na.rm=T)) %>% 
    mutate(ndc1_absolute_average = mean(c(
      ndc1_absolute_uncond, ndc1_absolute_cond), na.rm=T)) %>% 
    # Re-scale the targets so that higher values are cuts, zero is constant
    mutate_at(vars(ndc1_percentage_uncond, ndc1_percentage_cond, ndc1_percentage_average),
              ~(1-.)*100) %>% 
    # Get GHG growth rates
    mutate(ghg_growth_to_2015 = (ghg_cw_chosen_2015/ghg_cw_chosen_2000)^(1/15),
           ghg_growth_to_ndc1 = (ndc1_absolute_average/ghg_cw_chosen_2015)^(1/15)) %>% 
    mutate(growth_pp_pre = (ghg_growth_to_2015 - 1)*100,
           growth_pp_post_ndc1 = (ghg_growth_to_ndc1 - 1)*100,
           ndc1_delta_growth = growth_pp_post_ndc1 - growth_pp_pre) %>% 
    # Drop repeated observations
    select(-year, -ghg_primap, -ghg_cw_ex, -ghg_cw_lucf, -ghg_cw_chosen) %>% 
    distinct() %>% 
    filter(iso3c %in% core_iso3c) %>% 
    # Keep only the core variables
    select(iso3c, country, ghg_cw_chosen_2010, ghg_cw_chosen_2015,
           ndc1_absolute_average, ndc1_absolute_uncond, ndc1_absolute_cond,
           ndc1_percentage_uncond, ndc1_percentage_cond, ndc1_percentage_average,
           ndc1_delta_growth, growth_pp_post_ndc1, growth_pp_pre) %>%
    mutate_at(vars(ghg_cw_chosen_2010:growth_pp_pre), ~replace(., is.na(.), NA)) %>% 
    as.data.frame()
  
  # Merge 2021 targets with GHG time series
  targets2021 <- full_join(primap, cw, by = c("iso3c", "year")) %>% 
    full_join(., ndc2021 , by = "iso3c", relationship = "many-to-many") %>% 
    arrange(iso3c, year) %>% 
    # Get emissions in base year and GHG 2010 as own variables
    mutate(ndc2_baseyear_emissions = ifelse((cw_ndc2_ghgtarget_type == "Base Year" &
                                               year == cw_ndc2_ghgtarget_referenceyear),
                                            ghg_cw_chosen, NA),
           ghg_cw_chosen_2000 = ifelse(year == 2000, ghg_cw_chosen, NA),
           ghg_cw_chosen_2010 = ifelse(year == 2010, ghg_cw_chosen, NA),
           ghg_cw_chosen_2015 = ifelse(year == 2015, ghg_cw_chosen, NA)) %>% 
    group_by(iso3c) %>% 
    mutate(ndc2_baseyear_emissions = my.max(ndc2_baseyear_emissions),
           ghg_cw_chosen_2000 = my.max(ghg_cw_chosen_2000),
           ghg_cw_chosen_2010 = my.max(ghg_cw_chosen_2010),
           ghg_cw_chosen_2015 = my.max(ghg_cw_chosen_2015)) %>% 
    ungroup() %>% 
    # Get absolute emissions levels in targets
    ungroup() %>% 
    mutate(ndc2_absolute_uncond = 
             ifelse(cw_ndc2_ghgtarget_type == "Fixed Level", cw_ndc2_ghgtarget_fixedlevel_goal, # For fixed level targets
                    ifelse(cw_ndc2_ghgtarget_type == "Scenario", # For scenario
                           cw_ndc2_ghgtarget_scenario_referencelevel*(1-cw_ndc2_ghgtarget_scenario_unconditional_goal/100),
                           ifelse(cw_ndc2_ghgtarget_type == "Base Year", # For base year
                                  ndc2_baseyear_emissions*(1-cw_ndc2_ghgtarget_baseyear_goal/100), NA)))) %>% 
    mutate(ndc2_absolute_cond = 
             ifelse(cw_ndc2_ghgtarget_type == "Fixed Level", cw_ndc2_ghgtarget_fixedlevel_conditional,
                    ifelse(cw_ndc2_ghgtarget_type == "Scenario",
                           cw_ndc2_ghgtarget_scenario_referencelevel*(1-cw_ndc2_ghgtarget_scenario_goal_conditional/100),
                           ifelse(cw_ndc2_ghgtarget_type == "Base Year",
                                  ndc2_baseyear_emissions*(1-cw_ndc2_ghgtarget_baseyear_goal_conditional/100), NA)))) %>% 
    # Get emissions levels in targets as percentage of 2010 emissions (same benchmark as PRIMAP)
    mutate(ndc2_percentage_uncond = ndc2_absolute_uncond/ghg_cw_chosen_2010,
           ndc2_percentage_cond = ndc2_absolute_cond/ghg_cw_chosen_2010) %>% 
    # Get an average across both NDC targets (unconditional and conditional)
    rowwise() %>% 
    mutate(ndc2_percentage_average = mean(c(
      ndc2_percentage_uncond, ndc2_percentage_cond), na.rm=T)) %>% 
    mutate(ndc2_absolute_average = mean(c(
      ndc2_absolute_uncond, ndc2_absolute_cond), na.rm=T)) %>% 
    # Re-scale the targets so that higher values are cuts, zero is constant
    mutate_at(vars(ndc2_percentage_uncond, ndc2_percentage_cond, ndc2_percentage_average),
              ~(1-.)*100) %>% 
    ungroup() %>% 
    # Get GHG growth rates
    mutate(ghg_growth_to_2015 = (ghg_cw_chosen_2015/ghg_cw_chosen_2000)^(1/15),
           ghg_growth_to_ndc2 = (ndc2_absolute_average/ghg_cw_chosen_2015)^(1/15)) %>% 
    mutate(growth_pp_pre = (ghg_growth_to_2015 - 1)*100,
           growth_pp_post_ndc2 = (ghg_growth_to_ndc2 - 1)*100,
           ndc2_delta_growth = growth_pp_post_ndc2 - growth_pp_pre) %>% 
    # Drop repeated observations
    select(-year, -ghg_primap, -ghg_cw_ex, -ghg_cw_lucf, -ghg_cw_chosen) %>% 
    distinct() %>% 
    # Keep only "final" NDC 
    group_by(iso3c) %>% 
    mutate(entry = row_number(),
           entry_max = my.max(entry)) %>% 
    ungroup() %>% 
    filter(entry==entry_max) %>% 
    filter(iso3c %in% core_iso3c) %>%
    # Keep only the core variables
    select(iso3c, country, ndc2_date = date, ghg_cw_chosen_2010, 
           ndc2_absolute_average, ndc2_absolute_uncond, ndc2_absolute_cond,
           ndc2_percentage_uncond, ndc2_percentage_cond, ndc2_percentage_average,
           ndc2_delta_growth, growth_pp_post_ndc2) %>%
    mutate_at(vars(ghg_cw_chosen_2010:growth_pp_post_ndc2), ~replace(., is.na(.), NA)) %>% 
    as.data.frame()
}


#### Calculating average growth rates in countries with BAU targets ####


growth_df <- full_join(
  ndc2015 %>% select(iso3c, cw_ndc1_ghgtarget_scenario_referencelevel),
  ndc2021 %>% select(iso3c, cw_ndc2_ghgtarget_scenario_referencelevel),
  by = 'iso3c') %>% 
  full_join(., targets2015 %>% 
              select(iso3c, ghg_cw_chosen_2015, growth_pp_pre, growth_pp_post_ndc1),
            by = 'iso3c') %>% 
  full_join(., targets2021 %>% 
              select(iso3c, growth_pp_post_ndc2),
            by = 'iso3c') %>% 
  full_join(., data_for_models %>% select(iso3c, region, income),
            by = 'iso3c') %>% 
  mutate(growth_pp_to_ndc1_reference = 
           (((cw_ndc1_ghgtarget_scenario_referencelevel/ghg_cw_chosen_2015)^(1/15)-1)*100),
         growth_pp_to_ndc2_reference = 
           (((cw_ndc2_ghgtarget_scenario_referencelevel/ghg_cw_chosen_2015)^(1/15)-1)*100)) %>% 
  select(iso3c, region, income,
         growth_pp_pre, 
         growth_pp_to_ndc1_reference, growth_pp_to_ndc2_reference,
         growth_pp_post_ndc1, growth_pp_post_ndc2) %>% 
  filter(income != "High income") %>% 
  mutate(ratchet = factor(ifelse(growth_pp_to_ndc2_reference < growth_pp_to_ndc1_reference, 1, 0))) %>% 
  filter(!is.na(ratchet)) %>% 
  # Remove wild outliers
  filter(growth_pp_pre < 12, growth_pp_pre > -5) %>%
  filter(growth_pp_post_ndc1 < 20, growth_pp_post_ndc1 > -10) %>%
  filter(growth_pp_post_ndc2 < 25, growth_pp_post_ndc2 > -20)

## Calculate average growth rates
growth_df %>% 
  mutate(mean_pre = mean(growth_pp_pre)) %>% 
  group_by(ratchet) %>% 
  reframe(mean_pre, 
          mean_delta_reference_ndc1 = mean(growth_pp_to_ndc1_reference),
          mean_delta_reference_ndc2 = mean(growth_pp_to_ndc2_reference)) %>% 
  distinct() 
#' where "mean_pre" = mean growth rate for all countries in the pre-Paris period
#' where "mean_delta_reference_ndc1" = mean growth rate to NDC1's reference/BAU level
#' where "mean_delta_reference_ndc2" = mean growth rate to NDC2's reference/BAU level


